#include "Stdafx.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include "AEPD_Dist.h"
//#include "DataTypesII.h"
#include "specialfunctions.h"




using namespace std;
using namespace alglib;



#define MIN_NEG 0.0000001
#define EDGE 5

/* overload absolute function to take doubles */
double abs (double number){
	if (number >= 0){
		return number;
	}
	else return -1*number;
}

void AEPD_Fit(double * x, int N, double * params){

	// Take the array in x[0] - x[N-1] and put that into the vector data.
	vector<double> data(x, x+N);

	/*for (int i = 0; i < N; i++){
//		cout << "X[" << i << "] = " << x[i]<< endl;
		cout << "Before Sort data[" << i << "] = "<< data[i]<<endl;
	}*/


	/* Sorts the data: can be done in InArray.updateArray */
	double tmp;
	for(int i = 0;i < N-1 ;i++){
		for (int j = i + 1;j < N;j++){
			if (data[i] > data[j]){
				tmp  = data[i];
				data[i] = data[j];
				data[j] = tmp;
			}
		}
	}


	///sort(data.begin(),data.end());

	//for (int i = 0;i< N; i++){
	//	cout << "data sorted element " << i << " = " << data[i] << endl;
	//}

	//for (int i = 0;i< 4; i++){
	//	cout << "params " << i << " = " << params[i] << endl;
	//}
	double INCREMENT_ALPHA = 0.1;
	int sizeALPHA = 80;
	int NUM_STEPS = 500;

	double inc_A = INCREMENT_ALPHA;//////Iteration Increment of Alpha Parameter
	double start_A = 0.7;//////// Value to Start Alpha Iteration
	int sizeA = sizeALPHA;//// Number of iterations to perform

	double ** Kappa = new double*[sizeA];
	double ** Sigma = new double*[sizeA];
	double ** H = new double*[sizeA];
	for (int i = 0 ; i< sizeA ; i++){
		Kappa[i] = new double[N];
		Sigma[i] = new double[N];
		H[i] = new double[N];
	}

	double * X_minus = new double[sizeA];
	double * X_plus = new double[sizeA];


	int idx_k = 0;
	int idx_q = 0;
	double min_H = 10000000000;

	//int stepX = round(N/NUM_STEPS);
	int stepX = N/NUM_STEPS;
	int X_right = 0;
	int X_left = 0;
	if (stepX <1)
		stepX=1;

	for(int k  = EDGE -1 ; k<= N- EDGE -1 ; k=k+ stepX){

		double p1 = data[k];

		X_right = N- k;
		X_left = k ;

		double * x_minus_tmp = new double[X_left];
		double * x_plus_tmp = new double[X_right];

		/* Remember that data[] is already sorted.
		 * Take k as the current index (p1 = data[k]),
		 * x_minus_tmp[0 -> k] becomes absolute value (data - p1)
		 * x_plus_tmp[0 -> N - k] is just data -p1 (because already sorted)
		 */

		for( int q = 0; q< N; q++){
			if (q < X_left)
				x_minus_tmp[q] = abs(data[q] - p1);
			else
				x_plus_tmp[q - X_left] = data[q] - p1;			
		}

		/* start_A = 0.7
		 * q = [0 -> sizeA]
		 * inc_A = incrementAlpha
		 * alpha = 0.7 + q * inc_A
		 * Basically alpha goes from the start value through sizeA increments of inc_A 		 *
		 */


		for ( int q = 0; q < sizeA; q++){
			X_minus[q] = 0;
			X_plus[q] = 0;

			double alpha = start_A + q * inc_A;
			/* Goes through the array, If this is matlab then the code would be:
			 * alpha[sizeA] = 0.7 + sizeA*inc_A
			 * data[N] = data vector coming in
			 * k = index value			 *
			 * x_minus_tmp [k];
			 * x_plus_tmp[N-k];
			 * x_minus_tmp = abs( data[0-k] - data[k] );
			 * x_plus_tmp = data[k-N] - data[k];
			 * X_minus[q] = sum (x_minus_tmp.^alpha[q])/N;
			 * X_plus[q] = sum (x_plus_tmp.^alpha[q])/N;
			 *
			 * So the outer loop pins k, the inner loop pins q.
			 * and for each nested loop we have X_minus[k, q] and X_plus [k,q]
			 *
			 */

			for (int t = 0; t < N; t++){
				if (t < X_left)
					X_minus[q] += pow(x_minus_tmp[t], alpha) / N;
				else 
					X_plus[q] += pow(x_plus_tmp[t - X_left], alpha)/ N;
			}
		}


		//Produce potential parameter sets as a function of  Alpha, Sigma, Kappa which will be used for maximization
		for (int q = 0; q < sizeA; q++){

			double alpha = start_A + q * inc_A;
			Kappa[q][k] = pow( X_minus[q] / X_plus[q], 1 / (2*(alpha +1)));

			Sigma[q][k] = pow((alpha*  pow( X_minus[q]*X_plus[q], alpha/ (2*(alpha+1)))\
				* ( pow(X_plus[q], 1/(alpha+1)) + pow(X_minus[q], 1/(alpha+1)))), 1/alpha); 

			double sgngam = 0;
			H[q][k] = log(Sigma[q][k]) + lngamma(1/alpha,sgngam) \
				+ log(Kappa[q][k] + 1/Kappa[q][k]) + 1/alpha - log(alpha);

			//cout << "Kappa : " << Kappa[q* N + k] << "\t" << "Sigma : " << Sigma[q* N + k] << "\t"<< "H :" << H[q* N + k] << endl;
			// keep record of the minimum entropy index
			if(q < sizeA -2 && q >=1 && H[q][k] < min_H ){
				min_H = H[q][k];
				idx_k = k;
				idx_q = q;
			}
		}

		delete [] x_minus_tmp;
		delete [] x_plus_tmp;
	}

	params[0] = data[idx_k]; // theta
	params[1] = Sigma[idx_q][idx_k]; // sigma
	params[2] = Kappa[idx_q][idx_k]; // kappa
	params[3] = start_A + idx_q * inc_A; // alpha

	for (int i = 0 ; i< sizeA ; i++){
		delete [] Kappa[i];
		delete [] Sigma[i];
		delete [] H[i];

	}

	delete [] Kappa;
	delete [] Sigma;
	delete [] H;
	delete [] X_minus;
	delete [] X_plus;
	cout << "Got here with no problem "<< endl;
}


double AEPD_PDF(double * x, int N, double theta, double sigma, double kappa, double alpha){
	double g = 0;
	double D = 0;


	for (int i = 0;i< N ; i++)
	{
		if( x[i] >= theta)
			g = alpha/(sigma*gammafunction(1/alpha))* kappa/(1+pow(kappa,2)) \
			* exp( - pow(kappa, alpha)/pow(sigma,alpha)*pow( x[i]- theta, alpha));

		else
			g = alpha/(sigma*gammafunction(1/alpha))* kappa/(1+pow(kappa,2)) \
			* exp( -  1/( pow(kappa, alpha) *pow(sigma,alpha)) *pow( abs(x[i]- theta), alpha));
		if ( g<= 0)
			g = MIN_NEG;

		D =D+g;

	}

	return g/D;
}



// incompletegammac: "upper"
// incompletegamma: "lower"
double AEPD_CDF(double x, double theta,double sigma, double kappa, double alpha){


	double g = 0;
	double u = x - theta;

	if( u < 0){
		g = pow(kappa, 2) / ( 1+ pow(kappa,2)) * incompletegammac(1/alpha, pow(abs(u)/(sigma*kappa),alpha));
	}
	else
		g = 1 - 1/( 1+ pow(kappa,2))* incompletegammac(1/alpha, pow(u,alpha) * pow(kappa, alpha) / pow(sigma, alpha));

	return g;
}



// notice we use invincompletegammac here
// while MATLAB use the lower tail
double AEPD_Inv(double Percentile, double theta,double sigma,  double kappa, double alpha){


	double x = 0;
	if(Percentile < pow(kappa,2)/(1+pow(kappa,2)))
		x = theta-(1-exp( -sigma*kappa\
		*pow(invincompletegammac(1/alpha, 1/(pow(kappa,2)/(1+pow(kappa,2))) *Percentile), 1/alpha))) ;
	else
		x = theta+(exp( sigma/kappa*(pow (invincompletegammac( 1/alpha,(1-Percentile)*(1+pow(kappa,2))),1/alpha)))-1);

	return x;
}

void AEPD_Stat(double theta,double sigma,  double kappa, double alpha, double  *ret){
	ret[0] = (theta+sigma*(1/kappa-kappa)*gammafunction(2/alpha)/gammafunction(1/alpha));
	ret[1] = pow(sigma,2) * gammafunction(3/alpha)/gammafunction(1/alpha)* (1+pow(kappa,6))/(pow(kappa,2)+pow(kappa,4))\
		- pow(sigma,2) * pow(gammafunction(2/alpha),2)/ pow(gammafunction(1/alpha),2)* pow(1- pow(kappa,2),2)/ pow(kappa,2);	
}

void fast_AEPD_Fit(double * x, int N, double * params, struct DeviationStruct *D){}
/*
void fast_AEPD_Fit(double * x, int N, double * params, struct DeviationStruct *D){

	vector<double> data(x, x+N);
	sort(data.begin(),data.end());



	double inc_A = INCREMENT_ALPHA;//////Iteration Increment of Alpha Parameter
	double start_A = 0.6;//////// Value to Start Alpha Iteration
	int sizeA = sizeALPHA;//// Number of iterations to perform

	double Kappa[sizeALPHA][NUM_FITS*WINDOW_SIZE];
	double Sigma[sizeALPHA][NUM_FITS*WINDOW_SIZE];
	double  H[sizeALPHA][NUM_FITS*WINDOW_SIZE];
	////for (int i = 0 ; i< sizeA ; i++){
		///Kappa[i] = new double[N];
		///Sigma[i] = new double[N];
		///H[i] = new double[N];
	///}

	double  X_minus[sizeALPHA];
	double  X_plus[sizeALPHA];


	int idx_k = 0;
	int idx_q = 0;
	double min_H = 10000000000;

	int stepX = round(N/NUM_STEPS);
	int X_right = 0;
	int X_left = 0;
	if (stepX <=1)
		stepX=1;
	int x_minus_tmp[NUM_FITS*WINDOW_SIZE];
	int x_plus_tmp[NUM_FITS*WINDOW_SIZE];

	double p1;

	for(int k  = EDGE -1 ; k<= N- EDGE -1 ; k=k+ stepX){

		p1 = data[k];

		X_right = N- k-1;
		X_left = k+1 ;


		for( int q = 0; q< N; q++){
			if(q<=k){
				x_minus_tmp[q] = round(FACTOR*abs(data[q] - p1));
				x_plus_tmp[q]=0;
			}
			else if(q>k){
				x_plus_tmp[q] = round(FACTOR*abs(data[q] - p1));	
				x_minus_tmp[q]=0;
			}
		}


		for ( int q = 0; q < sizeA; q++){
			X_minus[q] = 0;
			X_plus[q] = 0;

			for( int t = 0; t< N; t++){
				if(x_minus_tmp[t]>0){
					///cout<<"x_minus:\t"<< D[x_minus_tmp[t]][q]<<endl;
					X_minus[q] += D[x_minus_tmp[t]].shapeIDX[q] / N;
				}
				if(x_plus_tmp[t]>0){	
					
					X_plus[q] += D[x_plus_tmp[t]].shapeIDX[q]/ N;
				}
			}
		}


		//Produce potential parameter sets as a function of  Alpha, Sigma, Kappa which will be used for maximization
		for (int q = 0; q < sizeA; q++){

			double alpha = start_A + q * inc_A;
			Kappa[q][k] = pow( X_minus[q] / X_plus[q], 1 / (2*(alpha +1)));

			Sigma[q][k] = pow((alpha*  pow( X_minus[q]*X_plus[q], alpha/ (2*(alpha+1)))\
				* ( pow(X_plus[q], 1/(alpha+1)) + pow(X_minus[q], 1/(alpha+1)))), 1/alpha); 

			double sgngam = 0;
			H[q][k] = log(Sigma[q][k]) + lngamma(1/alpha,sgngam) \
				+ log(Kappa[q][k] + 1/Kappa[q][k]) + 1/alpha - log(alpha);

			//cout << "Kappa : " << Kappa[q* N + k] << "\t" << "Sigma : " << Sigma[q* N + k] << "\t"<< "H :" << H[q* N + k] << endl;
			// keep record of the minimum entropy index
			if(q < sizeA -1 && q >=1 && H[q][k] < min_H ){
				min_H = H[q][k];
				idx_k = k;
				idx_q = q;
			}
		}

		///delete [] x_minus_tmp;
		///delete [] x_plus_tmp;
	}

	params[0] = data[idx_k]; // theta
	params[1] = Sigma[idx_q][idx_k]; // sigma
	params[2] = Kappa[idx_q][idx_k]; // kappa
	params[3] = start_A + idx_q * inc_A; // alpha

	/*for (int i = 0 ; i< sizeA ; i++){
		delete [] Kappa[i];
		delete [] Sigma[i];
		delete [] H[i];

	}

	delete [] Kappa;
	delete [] Sigma;
	delete [] H;
	delete [] X_minus;
	delete [] X_plus;

}
*/
